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Monte-Carlo study of the semimetal-insulator phase transition in monolayer graphene 
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We report on the results of the first-principle numerical study of spontaneous breaking of chiral 
(sublattice) symmetry in suspended monolayer graphene due to electrostatic interaction, which 
takes into account the screening of Coulomb potential by electrons on cr-orbitals. In contrast to 
the results of previous numerical simulations with unscreened potential, we find that suspended 
graphene is in the conducting phase with unbroken chiral symmetry. This finding is in agreement 
with recent experimental results by the Manchester group. Further, by artificially increasing the 
interaction strength we demonstrate that suspended graphene is quite close to the phase transition 
associated with spontaneous chiral symmetry breaking, which suggests that fluctuations of chirality 
and nonperturbative effects might still be quite important. 

PACS numbers: 73.22.Pr, 71.30.+h, 05.10.Ln 



In recent years significant effort has been invested into 
numerical studies of the electronic transport properties of 
ideal monolayer graphene [ ]. Since the electromagnetic 
coupling constant in graphene is effectively enhanced by 
the factor c/vp ~ 300, where c is the velocity of light 
and vf is Fermi velocity, charge carriers turn out to be 
strongly coupled, and various non-perturbative phenom- 
ena such as spontaneous breaking of chiral (sublattice) 
symmetry can emerge. In this situation first-principle 
calculations can only be performed numerically. 

In the seminal works [ ] it has been realized that 
graphene at neutrality point can be efficiently simulated 
by the Hybrid Monte-Carlo method, which is commonly 
used in lattice Quantum Chromodynamics (QCD). In 
the more recent work [3] Hybrid Monte- Carlo method 
was applied to perform a direct simulation of the tight- 
binding model of monolayer graphene (the possibility of 
such simulations was also discussed in [ ]). In these sim- 
ulations only the nearest-neighbour hopping for the tt 
orbitals was considered, and inter-electron interactions 
were described by the Coulomb law (with some finite on- 
site interaction potential). So far all simulations, both 
with the low-energy effective theory and with the tight- 
binding model, have indicated that at the critical cou- 
pling constant ac ~ 1 there is a semimetal-insulator 
phase transition associated with the emergence of a mass 
gap in the quasiparticle spectrum due to spontaneous chi- 
ral symmetry breaking. According to these results sus- 
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pended graphene, for which the effective coupling con- 
stant is as = e'^/hvF ^ 300/137 ^ 2.2, should be deeply 
in the insulating gapped phase with broken chiral sym- 
metry (we note also that in this phase graphene is in fact 
an anti- ferromagnetic [5]). 

However, these findings are in clear contradiction with 
recent experimental studies of the Manchester group [6] , 
in which no indications of the existence of a mass gap 
in suspended monolayer graphene were found. Till now 
the origin of this discrepancy between experimental and 
numerical data was not clear. In this paper we demon- 
strate that if one takes into account the screening of the 
Coulomb potential due to electrons on <j-orbitals of car- 
bon, the interaction between electrons should be even 
stronger than in suspended graphene in order to trig- 
ger the semimetal-insulator phase transition. To this 
end we perform Hybrid Monte-Carlo simulations of the 
tight-binding model of monolayer graphene with the par- 
tially screened inter-electron interaction potential ob- 
tained in [ ] in the constrained random phase approxi- 
mation (cRPA). In the calculations of [ ] only the screen- 
ing due to (7-orbitals was taken into account, thus one 
can use it as an input to the tight-binding model of 
electrons on 7r-orbitals without any double-counting of 
screening terms. The observed shift of the phase transi- 
tion thus eliminates the controversy between experimen- 
tal and numerical results and opens up the possibility of 
much more realistic first-principle Monte- Carlo simula- 
tions of the electronic properties of graphene. We fur- 
ther demonstrate that a rather mild increase of interac- 
tion strength do leads to spontaneous chiral symmetry 
breaking. Due to such proximity of the transition point, 
nonperturbative effects can be quite important in sus- 
pended graphene. 



Since the screening of the Coulomb potential due to 
cr-orbitals is mostly important at small distances of the 
order of lattice spacing [7] , it seems that the position of 
the semimetal-insulator phase transition is highly sensi- 
tive to the form of the inter-electron interaction potential 
at short distances. We note that the high sensitivity of 
low-energy effective theory to ultraviolet regularization 
was also discovered in the work [ ], where fermionic prop- 
agators were found to be saturated by momenta of the 
order of inverse lattice spacing. 

The fact that in suspended monolayer graphene the ef- 
fective inter-electron interaction should be weaker than 
in the tight-binding model for the tt orbitals was also 
noted in [9] by fitting the numerical value of the renor- 
malized Fermi velocity vf (<^) to the experimental data 
of [6]. The corresponding value of a was estimated as 
a ~ 0.7. ..0.9, which is significantly smaller than a^. 
A recent semi- analytic study [ ] based on Schwinger- 
Dyson equations has also shown that the phase transition 
is shifted to larger couplings if one takes into account the 
renormalization of the Fermi velocity. 

The starting point of our simulations is the tight- 
binding Hamiltonian with the staggered potential m: 

Htb = -^ X] (ctlctx + blbx + h.cA + 

<x,y> 

+ \^±maj.aa; ± m6j.6a;. (1) 
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FIG. 1: A comparison of the partially screened Coulomb po- 
tential with the exact Coulomb potential and the potential 
obtained from non-compact gauge field on the hexagonal lat- 
tice [3]. 



responds to the Dirac mass. 

Next we introduce the interaction Hamiltonian with an 
inter-electron interaction potential Vxy'. 



He = - ^VxyqxQy, 



(2) 



where tz = 2.7 eV, the sum ^ is performed over all 

<x,y> 

pairs of nearest-neighbour sites of the graphene hexag- 
onal lattice (we impose periodic spatial boundary con- 
ditions as in [3]) and d^ , a and b^ , b are the cre- 
ation/annihilation operators for particles and holes, re- 
spectively. The latter are related to creation/annihilation 
operators cj, g, Cx^s ^oi electrons with spin s =t,i as 

o^x = Cx^^j bx = ±c^ I , where we take the plus sign for 
X belonging to one of the simple sublattices of graphene 
hexagonal lattice and the minus sign - for another simple 
sublattice [3, ]. The whole Hilbert space of the tight- 
binding model can be constructed by the action of the 
creation operators aj., b^ on the ground state |0) fixed 
by the conditions a^^ |0) = 0, 6a; |0) = 0. In this ground 
state each lattice site is occupied by one electron with 
spin down. Of course, in Monte-Carlo simulations we 
sum over all possible states of the system, so this choice 
of the ground state is only motivated by calculational 
convenience. 

The staggered potential is equal to -\-m for the sites of 
one simple sublattice and —m for sites of another sim- 
ple sublattice. Its role is twofold: first, it regularizes the 
inverse of the fermionic kinetic operator in the Hybrid 
Monte-Carlo algorithm [3, 4]. Second, the staggered po- 
tential explicitly breaks the chiral (sublattice) symmetry 
and thus serves as a seed for spontaneous chiral symme- 
try breaking, which would otherwise be impossible in a 
finite volume. In the low-energy effective theory m cor- 



where qx = dl^dx 
at lattice site x. 



blbx is the operator of electric charge 



For the on-site interaction potential Vxx = ^oo and 
the potentials between nearest (l^i), next-to- nearest 
(F02) and next-to- next-to- nearest-neighbouring lattice 
sites (F03) we use the values calculated in [7] (see ta- 
ble I, 3d column). The resulting shape of the potential 
is illustrated on Fig. 1. At larger distances we use the 
Coulomb potential V (r) = 1/ (car). The form of the po- 
tential is additionally corrected to account for periodic 
boundary conditions. The factor e^r ~ 1.41 describes the 
contribution of electrons on a orbitals to the effective di- 
electric permittivity of graphene monolayer at large dis- 
tances and is obtained by equating Vqs to the Coulomb 
potential at r = ros = 0.284 nm: F03 = 1/ (ecr^os)- Phys- 
ically this means that we assume that all the charges 
which screen the potential of a test charge are localized 
within the radius ro3. It is important to stress that this 
large-distance correction of the potential by a factor l/ccr 
alone is insufficient to prevent the semimetal-insulator 
phase transition in suspended graphene. Indeed, since 
for the unscreened Coulomb potential the corresponding 
critical value of the coupling constant ac ~ 1 [ , ] is 
more than two times smaller than the effective coupling 
constant ag ~ 2.2 in suspended graphene, the coefficient 
Co- should be at least larger than 2 in order to shift the 
phase transition to a > as . 

We proceed by making the standard Suzuki- Trotter 



decomposition of the partition function: 



Nt 



= Tr Te"^" 



S^-HcS-HtbS 



...)+^(^')^ (3) 



where /3 = {kT)~ is the inverse temperature and S = 
P/Nt with TVt > 1. The factors in the last hue of (3) 
are now interleaved with decompositions of the identity 
operator over Grassmann coherent states: 

1=1 dipdrjdipdf] e ^ ^ \ip ^ r]) [ip ^ r]\ ^ 



|V^,7/) 



|0). (4) 



The matrix elements (ip^rjle ^^*^ \ip' ^r]') can be now eas- 
ily calculated using the identity 



/ _i ^x "^y y 



m=exJj2^^{e^)^J'\. (5) 



.x,y 



In order to find the matrix elements of the exponent 
of the interaction Hamiltonian He we perform the 
Hubbard- Stratonovich transformation [4]: 

\ x,y / 

^ / V(f^exp i-^^^^xV-y^ify -iS'^if^q^j, (6) 

\ x,y X / 

where V^^ is the matrix inverse of the potential Vxy: 
^ ^xz^^zy = ^xy After that we again apply the formula 

z 

(5) to the last line of (6) and finally arrive at the fol- 
lowing functional integral representation of the partition 
function: 



Tre 



-(5H 



/p^.,.p^.,.p,.,.p^.,.p,.,. 



,(7) 



where S [^x,n] = f XI '^x,nV~^Lpy^n is the action of the 

x^y^n 

Hubbard field (px,n and n = ... 2Nt — 1 enumerates 
the factors in the last line of (3). The fermionic part of 
the action is written as follows: 



/ ^ Y X ,n -^^-^ X ,y ,n ,n' 'y y ,n' 



Nt-1 

= E 



/e=0 



y^ 'lpx,2k {i^x,2k - V^x,2/c+l) 
_ X 

<x,y> 

X 

+ (5 ^ ±m-0a;,2fc^a;,2fe+l 



(8) 



In this expression the Grassmann variables ^l^x,2k and 
i^x,2k+i label the fermionic coherent states inserted be- 
tween the factors e"^*^^, e"^^^ and e"^^^, e"^*^^ in 
(3), respectively. It can be shown that such a "double- 
layer" structure of the action leads to discretization er- 
rors of the order of (5, in contrast to simpler fermionic 
action constructed in [4], for which discretization errors 
scale as V^. We also impose anti-periodic boundary con- 
ditions in time direction on fermionic variables ipx.n-, Vx.n 
which amounts to replacing (j)x,Nt-i ~^ ^x,Nt-i + 'tt in 
the fermionic action (8). 

Now the Grassmann variables in (7) can be integrated 
out, which yields the following representation of the par- 
tition function: 

T,e-P"^ /p^^_„e-Sbx,n]|det(M[^,,„])|2. (9) 



The manifest positivity of the integration weight in (9) 
is due to the symmetry between particles and holes for 
graphene at neutrality point. For example, at finite 
chemical potential the two fermionic determinants ap- 
pearing in (9) after integration over il)x,n and rjx^n iii 
(7) would no longer be complex conjugate, which would 
make Monte- Carlo simulations much more difficult due 
to the fermionic sign problem. For our choice of the inter- 
electron interaction potential, the action of the Hubbard 
field 5* [<^x,n] is also a positive definite quadratic form. 
Thus we can generate the configurations of (px,n by a 
Monte- Carlo method and calculate physical observables 
as averages over the generated configurations. Here we 
follow [3, 4] and use the Hybrid Monte-Carlo method 
with the <l>-algorithm. Inversion of fermionic operator M 
is the most difficult part of this algorithm, so CPUs were 
used to accelerate the calculations. 

In order to detect the chiral symmetry breaking, we 
calculate the chiral condensate, which is the difference of 
particle numbers between the two simple sublattices A 
and B: 

( An) = 1( ^(ata, + b%) - E(4«=. + b%) )(10) 



xeA 



xeB 



where N is the overall number of sites of one sublattice 
of hexagonal lattice. In terms of the fermionic operator 
Mx,y,n,n' this expectation value reads: 



(An) 



NNt 



2Nt-l 

E 

n=0 xeA 



xeB 



where the average is now taken over configurations of the 
Hubbard field with the weight (9). 

Our simulations were performed on the lattice with 
spatial size 18 x 18 and Nt = 20, J = 0.1eV"\ which 
corresponds to the temperature T = 0.5 eV = 5.8 • 10^ K. 
This temperature is considerably higher than in real ex- 
periments, however, in our simulations it is the temper- 
ature of the electron gas only. We do not consider ther- 
mal fiuctuations of the crystalline lattice, thus phonon 
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FIG. 2: The dependence of the chiral condensate (11) on e 
and on vn (in the inset) for the 18 x 18 lattice with Nt — 20 
and 8 — 0.1 eV~^. For e = 1.0 we show the results obtained 
on the 24 x 24 lattice with Nt = 40, ^ = 0.05 eV"\ 



temperature is formally zero. We rely here on the re- 
sults of [3], which indicate that as long as the electron 
temperature is much smaller than the hopping parame- 
ter K in (1), it does not significantly affect the insulator- 
semimetal phase transition. To study the behavior of the 
condensate (10) in the limit tti ^ 0, we perform simu- 
lations at five different values of the staggered potential: 
vn = 0.1, 0.2, 0.3, 0.4, 0.5 eV. The interaction strength 
is controlled by additionally rescaling the potential by 
some factor e: V^y -^ Vxy/e. 

The coefficient e can be thought of as the dielectric 
permittivity of the medium surrounding the graphene 
monolayer. However, to make this interpretation physi- 
cally consistent one should also redo the calculations of 
[7] taking into account this additional screening. In our 
case e has no direct physical interpretation and is only 
used to characterize the proximity of suspended graphene 
(which corresponds to e = 1) to the phase transition. For 
each set of lattice parameters we have generated 100 sta- 
tistically independent configurations of the field (px,n- 

The dependence of the chiral condensate (11) on e for 
6 < 1 is illustrated on Fig. 2. To obtain the plotted 
values of An, we have fitted the mass dependence of the 
condensate An (m) by a quadratic function of m and 
used this fit to extrapolate An (m) to m = 0. These fits 
are shown on Fig. 2 in the inset. One can see that the 
extrapolated value An {m -^ 0) for suspended graphene 
(e = 1) is equal to zero within error range, which indi- 



cates the absence of chiral symmetry breaking. We have 
also checked this result on the larger (24 x 24, Nf = 20, 
S = 0.1 eV"^) and finer (24 x 24, Nt = 40, ^ = 0.05 eV"^) 
lattices and on the larger set of 250 configurations of 
(px,n- All our measurements confirm that after extrapo- 
lation to m = the chiral condensate is equal to zero for 
suspended graphene. 

Only at e < ec ~ 0.7 the extrapolation to m ^ yields 
nonzero chiral condensate, which suggests that the state 
with broken chiral symmetry is favoured, and sponta- 
neous chiral symmetry breaking is likely in the infinite 
volume limit. The fact that the critical value Cc ^ 0.7 is 
quite close to one suggests that while suspended graphene 
is still in the conducting phase with unbroken chiral sym- 
metry, the proximity of the phase transition can still 
manifest itself in large fluctuations of order parameter 
(chiral condensate) and in other non-perturbative phe- 
nomena. 

We conclude that the screening of the Coulomb po- 
tential by electrons on a-orbitals strongly influences 
the insulator-semimetal phase transition in monolayer 
graphene, so that the transition point is shifted into 
the region of parameter space in which the interaction 
strength is even stronger than in suspended graphene. 
We also note an intriguing possibility to effectively en- 
hance the inter-electron interactions by stretching the 
graphene layer [ ], which can be used to reach the tran- 
sition point in experiment. 
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